Load packages
library(Matrix)
library(data.table)
library(Seurat)
library(SeuratDisk)
library(dplyr)
library(Seurat)
dataset_name = "cell2location34"
sc_dataset_file = paste(dataset_name,"sc_dataset.rds",sep ="_")
rds_file= sc_dataset_file
sc_seu <- readRDS(rds_file)
rds_file='Visium_spatial.rds'
sp_seu <- readRDS(rds_file)
if(file.exists("old_version_ls.rds")){
rds_file= "old_version_ls.rds"
spotlight_ls <- readRDS(rds_file)
nmf_mod <- spotlight_ls[[1]]
decon_mtrx <- spotlight_ls[[2]]
}else{
sc_seu <- SCTransform(sc_seu, verbose = FALSE) %>% RunPCA(verbose = FALSE)
sp_seu <- SCTransform(sp_seu,assay = "Spatial", verbose = FALSE) %>% RunPCA(verbose = FALSE)
#### Extract the top marker genes from each cluster ####
Seurat::Idents(object = sc_seu) <- sc_seu@meta.data$Subset
cluster_markers_all <- Seurat::FindAllMarkers(object = sc_seu,
assay = "RNA",
slot = "data",
verbose = TRUE,
only.pos = TRUE,
logfc.threshold = 1,
min.pct = 0.9)
set.seed(123)
setwd("/data/home/lyang/Visium_spotlight")
start_time <- Sys.time()
spotlight_ls <- spotlight_deconvolution(se_sc = sc_seu,
counts_spatial = sp_seu@assays$Spatial@counts,
clust_vr = "Subset",
cluster_markers = cluster_markers_all,
cl_n = 100, # 100 by default
hvg = 5000,
ntop = NULL,
transf = "uv",
method = "nsNMF",
min_cont = 0.09)
saveRDS(spotlight_ls,"old_version_ls.rds")
end_time <- Sys.time()
end_time - start_time
}
unloadNamespace("SPOTlight")
library(SPOTlight, lib.loc="/data/home/bioinfo/R/R4.1.0/library_B3.13/old_versions/")
h <- NMF::coef(nmf_mod[[1]])
rownames(h) <- paste("Topic", 1:nrow(h), sep = "_")
topic_profile_plts <- SPOTlight::dot_plot_profiles_fun(
h = h,
train_cell_clust = nmf_mod[[2]])
topic_profile_plts[[2]] + ggplot2::theme(
axis.text.x = ggplot2::element_text(angle = 90),
axis.text = ggplot2::element_text(size = 12))

topic_profile_plts[[1]] + theme(axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 12))

# change the cell type names in sp_seu, replace _ with . to keep names coincide between sp_seu and decon_mtrx
decon_mtrx_sub <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]
decon_mtrx_sub[decon_mtrx_sub < 0.08] <- 0
decon_mtrx <- cbind(decon_mtrx_sub, "res_ss" = decon_mtrx[, "res_ss"])
rownames(decon_mtrx) <- colnames(sp_seu)
decon_df <- decon_mtrx %>%
data.frame() %>%
tibble::rownames_to_column("barcodes")
sp_seu@meta.data <- sp_seu@meta.data %>%
tibble::rownames_to_column("barcodes") %>%
dplyr::left_join(decon_df, by = "barcodes") %>%
tibble::column_to_rownames("barcodes")
# Seurat::SpatialFeaturePlot(
# object = sp_seu,
# features = c("B.Cycling", "B.GC.DZ"),
# alpha = c(0.1, 1))
cell_types_all <- colnames(decon_mtrx)[which(colnames(decon_mtrx) != "res_ss")]
getwd()
## [1] "/data/home/lyang/Visium_spotlight"
# install.packages("imager")
# SPOTlight::spatial_scatterpie(se_obj = sp_seu,
# cell_types_all = cell_types_all,
# img_path = "tissue_lowres_image.png",
# pie_scale = 0.4)
#
#
# SPOTlight::spatial_scatterpie(se_obj = sp_seu,
# cell_types_all = cell_types_all,
# img_path = "tissue_lowres_image.png",
# cell_types_interest = "B.naive",
# pie_scale = 0.4)
# "T_CD4+_naive","FDC","B_naive"
# This can be dynamically visualized with DT as shown below
unloadNamespace("SPOTlight")
library(SPOTlight)
ct <- colnames(decon_mtrx)
decon_mtrx[decon_mtrx < 0.1] <- 0
library(RColorBrewer)
if(dataset_name == "Tabula"){
n <- 29
}else if(dataset_name == "cell2location"){
n <- 44
}else if(dataset_name == "cell2location34"){
n <- 34
}
colrs <- brewer.pal.info[brewer.pal.info$colorblind == TRUE, ]
col_vec = unlist(mapply(brewer.pal, colrs$maxcolors, rownames(colrs)))
set.seed(123)
col <- sample(col_vec, n)
# show T cell and
# B cell zones and GCs with follicular dendritic cells (FDCs)
# FDC, T_CD4+_naive, B_naive
# Transcriptionally fine subtypes (B_GC_DZ, B_GC_LZ, B_GC_prePB and
# T_CD4+_TfH_GC); transcriptionally distinct subtypes (B_Cycling and FDC)
paletteMartin <- col
pal <- colorRampPalette(paletteMartin)(length(ct))
names(pal) <- ct
pal_back <- pal
plot_3_region <- function(pal)
{
for (char in names(pal)) {
print(char)
if(char %in% c("T.CD4..naive","FDC","B.naive") ){
if(char == "T.CD4..naive")
{
pal[char] = "#FFFF00"
} else if (char == "FDC")
{
pal[char] = "#add8e6"
} else if (char == "B.naive")
{
pal[char] = "#FF0000"
}
next
}
pal[char] = "#00008b"
}
return(pal)
}
packageVersion("SPOTlight")
## [1] '0.99.4'
pal = plot_3_region(pal_back)
## [1] "B.activated"
## [1] "B.Cycling"
## [1] "B.GC.DZ"
## [1] "B.GC.LZ"
## [1] "B.GC.prePB"
## [1] "B.IFN"
## [1] "B.mem"
## [1] "B.naive"
## [1] "B.plasma"
## [1] "B.preGC"
## [1] "DC.CCR7."
## [1] "DC.cDC1"
## [1] "DC.cDC2"
## [1] "DC.pDC"
## [1] "Endo"
## [1] "FDC"
## [1] "ILC"
## [1] "Macrophages.M1"
## [1] "Macrophages.M2"
## [1] "Mast"
## [1] "Monocytes"
## [1] "NK"
## [1] "NKT"
## [1] "T.CD4."
## [1] "T.CD4..naive"
## [1] "T.CD4..TfH"
## [1] "T.CD4..TfH.GC"
## [1] "T.CD8..CD161."
## [1] "T.CD8..cytotoxic"
## [1] "T.CD8..naive"
## [1] "T.TfR"
## [1] "T.TIM3."
## [1] "T.Treg"
## [1] "VSMC"
## [1] "res_ss"
library(ggplot2)
plotSpatialScatterpie(
x = sp_seu,
y = decon_mtrx,
cell_types = colnames(y),
img = FALSE,
scatterpie_alpha = 1,
pie_scale = 0.4) +
scale_fill_manual(
values = pal,
breaks = names(pal))

plot_n_cell_type <- function(pal,col_vec,n, cell_type_vector)
{
for (char in names(pal)) {
print(char)
if(char %in% cell_type_vector){
next
}
pal[char] = "#00008b"
}
return(pal)
}
cell_type_vector = c("B.GC.DZ", "B.GC.LZ", "B.GC.prePB","T.CD4..TfH.GC",
"B.Cycling", "FDC")
pal = plot_n_cell_type(pal_back,col_vec,6,cell_type_vector)
## [1] "B.activated"
## [1] "B.Cycling"
## [1] "B.GC.DZ"
## [1] "B.GC.LZ"
## [1] "B.GC.prePB"
## [1] "B.IFN"
## [1] "B.mem"
## [1] "B.naive"
## [1] "B.plasma"
## [1] "B.preGC"
## [1] "DC.CCR7."
## [1] "DC.cDC1"
## [1] "DC.cDC2"
## [1] "DC.pDC"
## [1] "Endo"
## [1] "FDC"
## [1] "ILC"
## [1] "Macrophages.M1"
## [1] "Macrophages.M2"
## [1] "Mast"
## [1] "Monocytes"
## [1] "NK"
## [1] "NKT"
## [1] "T.CD4."
## [1] "T.CD4..naive"
## [1] "T.CD4..TfH"
## [1] "T.CD4..TfH.GC"
## [1] "T.CD8..CD161."
## [1] "T.CD8..cytotoxic"
## [1] "T.CD8..naive"
## [1] "T.TfR"
## [1] "T.TIM3."
## [1] "T.Treg"
## [1] "VSMC"
## [1] "res_ss"
plotSpatialScatterpie(
x = sp_seu,
y = decon_mtrx,
cell_types = colnames(y),
img = FALSE,
scatterpie_alpha = 1,
pie_scale = 0.4) +
scale_fill_manual(
values = pal,
breaks = names(pal))

pal = pal_back
# type(pal)
plotSpatialScatterpie(
x = sp_seu,
y = decon_mtrx,
cell_types = colnames(y),
img = FALSE,
scatterpie_alpha = 1,
pie_scale = 0.4) +
scale_fill_manual(
values = pal,
breaks = names(pal))

rds_file='predictions.assay.rds'
predictions.assay <- readRDS(rds_file)
mat = decon_mtrx
spatial = sp_seu
predictions.assay@data <- t(mat)
meta.features <- as.data.frame(colnames(mat) )
rownames(meta.features) = meta.features[,1]
meta.features$`colnames(mat)` =NULL
predictions.assay@meta.features = meta.features
spatial[["predictions"]] <- predictions.assay
DefaultAssay(spatial) <- "predictions"
# table(sc_dataset@meta.data$Subset)
if(dataset_name == "Tabula"){
p = SpatialFeaturePlot(spatial, features = c("stromal cell", "t cell"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
# print(p)
# p = SpatialFeaturePlot(spatial, features = c("naive b cell", "b cell"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
# print(p)
}else if(dataset_name == "cell2location"){
SpatialFeaturePlot(spatial, features = c("T_Treg", "DC_cDC2"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
# SpatialFeaturePlot(spatial, features = c("DC_cDC1", "DC_cDC2"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)
}else if(dataset_name == "cell2location34"){
# detach("package:SpatialExperiment", unload=TRUE)
p = SpatialFeaturePlot(spatial, features = c("B.GC.LZ", "T.CD4..TfH.GC","B.GC.prePB","FDC"), pt.size.factor = 1.6, ncol = 4, crop = TRUE) + ggtitle("germinal center light zone")
print(p)
# + scale_fill_continuous(limits = c(0, 1))
p = SpatialFeaturePlot(spatial, features = c("B.Cycling", "B.GC.DZ"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)+ ggtitle("germinal center dark zone")
print(p)
SpatialFeaturePlot(spatial, features = c("B.naive", "B.preGC"), pt.size.factor = 1.6, ncol = 2, crop = TRUE) + ggtitle("B follicle + pre GC")
}



# table(sc_dataset@meta.data$Subset)